Twin analysis using TwinAnalysis

Ivan Voronin

12 января, 2021

Contents

Twin method: Fundamentals

ACE

ADE

Univariate ACE

Univariate ACE

Univariate ACE

Univariate ACE

Univariate ACE

TwinAnalysis

Univariate toy data

# Making univariate toy data
library(dplyr,
        quietly = TRUE)

mz_data <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0),
  Sigma = matrix(c(1.0, 0.8,
                   0.8, 1.0),
                 nrow = 2)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'X2'))

dz_data <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0),
  Sigma = matrix(c(1.0, 0.5,
                   0.5, 1.0),
                 nrow = 2)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'X2'))

data <- rbind(data.frame(zyg = 'MZ', mz_data),
              data.frame(zyg = 'DZ', dz_data)) %>%
  mutate(
    zyg = factor(zyg, levels = c('MZ', 'DZ'))
  )

head(data)
##   zyg          X1          X2
## 1  MZ  0.55107062 -0.54415443
## 2  MZ -1.77073118 -1.81837605
## 3  MZ  1.45332458  1.40457435
## 4  MZ -0.28017613  0.17413459
## 5  MZ -0.55666668 -1.33562604
## 6  MZ -0.03354105  0.01393325

Descriptives

# Installing the packages
# install.packages('devtools')
# library(devtools)
# install_github('ivanvoronin/mlth.data.frame')
# install_github('ivanvoronin/TwinAnalysis')

# Loading the packages
library(OpenMx)
library(TwinAnalysis)

# Descriptive statistics
get_univ_descriptives(data, zyg = 'zyg', vars = 'X')
## $X
##    Twin1-------------------- Twin2------------------- Cor--------------------------
##    n   M           SD        n   M           SD       r         lowerCI   upperCI  
## MZ 100  0.10347156 1.0247719 100  0.02369611 1.015420 0.8394995 0.7700587 0.8892814
## DZ 100 -0.01673161 0.9709349 100 -0.14933507 1.038595 0.5738636 0.4253948 0.6922529

Univariate ACE

# Univariate ACE model
model <- univariate_ace(data, zyg = 'zyg', ph = 'X')

model <- mxRun(model)

summary(model)
## Summary of ACE 
##  
## free parameters:
##            name matrix row col   Estimate  Std.Error A
## 1 ACE.mean[1,1]   mean   1   1 -0.0154077 0.06578862  
## 2    ACE.a[1,1]      a   X   X  0.7347096 0.08905736  
## 3    ACE.c[1,1]      c   X   X  0.5603042 0.12403442  
## 4    ACE.e[1,1]      e   X   X  0.4102589 0.02890131  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:              4                    396
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:               984.723
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/400
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:        192.723                992.723
## BIC:      -1113.411               1005.916
##       |  Sample-Size Adjusted
## AIC:                 992.9281
## BIC:                 993.2438
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:50:48 
## Wall clock time: 0.07298517 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Output - parameters

# Output tables from the univariate ACE
get_output_tables(model)
## $Variance_components
##        A(%)      C(%)      E(%)
## X 0.5281518 0.3071673 0.1646809
## 
## $Raw_variance
##      Total         A         C         E
## X 1.022051 0.5397982 0.3139407 0.1683124

Output - model fit

# Model fit
model <- ref_models(model)

fit_stats(model)
##   model      base ep minus2LL  df     AIC       BIC
## 1   ACE Saturated  4  984.723 396 192.723 -1113.411
##        CFI      TLI RMSEA   diffLL diffdf         p
## 1 1.000588 1.000196     0 5.905993      6 0.4338028

Constrained models

# Constrained models
# This is a list of MxModels
nested_models <- def_nested_models(
  model,
  AE = mxConstraint(c[1, 1] == 0),
  CE = mxConstraint(a[1, 1] == 0),
  E = list(mxConstraint(c[1, 1] == 0),
           mxConstraint(a[1, 1] == 0)),
  run = TRUE
)

fit_stats(model, nested_models = nested_models)
## Warning in computeFitStatistics(likelihood, DoF, chi,
## chiDoF, retval[["numObs"]], : Your model may be mis-
## specified (and fit worse than an independence model),
## or you may be using the wrong independence model, see ?
## mxRefModels
##   model      base ep  minus2LL  df      AIC        BIC
## 1   ACE Saturated  4  984.7230 396 192.7230 -1113.4107
## 2    AE       ACE  4  989.1442 397 195.1442 -1114.2878
## 3    CE       ACE  4 1007.3842 397 213.3842 -1096.0478
## 4     E       ACE  4 1144.5469 398 348.5469  -964.1834
##         CFI       TLI      RMSEA     diffLL diffdf
## 1 1.0005878 1.0001959 0.00000000   5.905993      6
## 2 0.9791975 0.9940564 0.04875019   4.421220      1
## 3 0.8651567 0.9614733 0.12411747  22.661211      1
## 4 0.0138350 0.7534587 0.31397644 159.823929      2
##              p
## 1 4.338028e-01
## 2 3.549466e-02
## 3 1.932294e-06
## 4 1.970946e-35

Bivariate toy data

# Making bivariate toy data
mz_data2 <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0, 0, 0),
  Sigma = matrix(
    c(1.0, 0.7, 0.6, 0.5,
      0.7, 1.0, 0.5, 0.8,
      0.6, 0.5, 1.0, 0.7,
      0.5, 0.8, 0.7, 1.0),
    nrow = 4)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'Y1', 'X2', 'Y2'))

dz_data2 <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0, 0, 0),
  Sigma = matrix(
    c(1.0, 0.7, 0.4, 0.2,
      0.7, 1.0, 0.2, 0.4,
      0.4, 0.2, 1.0, 0.7,
      0.2, 0.4, 0.7, 1.0),
    nrow = 4)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'Y1', 'X2', 'Y2'))

data2 <- rbind(data.frame(zyg = 'MZ', mz_data2),
               data.frame(zyg = 'DZ', dz_data2)) %>%
  mutate(
    zyg = factor(zyg, levels = c('MZ', 'DZ'))
  )

head(data2)
##   zyg         X1         Y1         X2          Y2
## 1  MZ  2.9737961  1.7046050  2.6551091  0.93346789
## 2  MZ -0.2515991 -0.1794637 -0.8486091 -0.65641665
## 3  MZ -1.2586533 -0.5918007 -1.1710607 -0.08288749
## 4  MZ  0.4271026  0.3153837  0.6676022  0.78442769
## 5  MZ  0.2249155  0.5210256 -0.4747569 -1.04644136
## 6  MZ  0.3755235  0.4408710 -0.4999070 -0.53263830

Descriptives

# Descriptive statistics
get_univ_descriptives(data2, 
                      zyg = 'zyg', 
                      vars = c('X', 'Y'))
## $X
##    Twin1------------------ Twin2------------------- Cor--------------------------
##    n   M          SD       n   M          SD        r         lowerCI   upperCI  
## MZ 100 0.06760087 1.146399 100 0.02576845 0.9835557 0.6359860 0.5023170 0.7399680
## DZ 100 0.01026312 1.000688 100 0.01358189 1.1121019 0.4952969 0.3310888 0.6303876
## 
## $Y
##    Twin1-------------------- Twin2-------------------- Cor--------------------------
##    n   M           SD        n   M           SD        r         lowerCI   upperCI  
## MZ 100 -0.01586990 0.9520076 100  0.02848424 0.9729878 0.7641308 0.6679673 0.8351952
## DZ 100 -0.01175412 1.0119057 100 -0.02396046 1.0687223 0.5158921 0.3554967 0.6467724

Twin correlations

# Cross-twin cross-trait correlations
get_cross_trait_cors(data2,
                     zyg = 'zyg',
                     vars = c('X', 'Y'))
## zyg: MZ
##           X1        Y1        X2        Y2
## X1 1.0000000 0.7778976 0.6359860 0.4868070
## Y1 0.7778976 1.0000000 0.5482844 0.7641308
## X2 0.6359860 0.5482844 1.0000000 0.7080632
## Y2 0.4868070 0.7641308 0.7080632 1.0000000
## -------------------------------------------- 
## zyg: DZ
##           X1        Y1        X2        Y2
## X1 1.0000000 0.7242958 0.4952969 0.3904858
## Y1 0.7242958 1.0000000 0.2686450 0.5158921
## X2 0.4952969 0.2686450 1.0000000 0.7006547
## Y2 0.3904858 0.5158921 0.7006547 1.0000000

Bivariate ACE

# Bivariate ACE model
model2 <- multivariate_ace(
  data = data2,
  zyg = 'zyg',
  ph = c('X', 'Y')
)

model2 <- mxRun(model2)

summary(model2)
## Summary of ACE 
##  
## free parameters:
##             name matrix row col     Estimate  Std.Error A
## 1  ACE.mean[1,1]   mean   1   1  0.028582736 0.06602367  
## 2  ACE.mean[1,2]   mean   1   2 -0.006697648 0.06398131  
## 3     ACE.a[1,1]      a   X   X  0.229633789 0.19036365  
## 4     ACE.a[2,1]      a   Y   X  0.487159593 0.14552821  
## 5     ACE.a[2,2]      a   Y   Y  0.761979977 0.09835996  
## 6     ACE.c[1,1]      c   X   X  0.525658510 0.07115624  
## 7     ACE.c[2,1]      c   Y   X  0.367411377 0.20938628  
## 8     ACE.c[2,2]      c   Y   Y  0.456055458 0.15872435  
## 9     ACE.e[1,1]      e   X   X  0.417842833 0.02948707  
## 10    ACE.e[2,1]      e   Y   X  0.493955852 0.05375753  
## 11    ACE.e[2,2]      e   Y   Y  0.467266109 0.03298510  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:             11                    789
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:              1793.504
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/800
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:       215.5039               1815.504
## BIC:     -2386.8685               1851.785
##       |  Sample-Size Adjusted
## AIC:                 1816.908
## BIC:                 1816.936
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:50:49 
## Wall clock time: 0.077667 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Output - parameters

# Output tables from the bivariate ACE
get_output_tables(model2)
## $Variance_components
##        A(%)      C(%)      E(%)
## X 0.2589905 0.3672559 0.3737536
## Y 0.5766131 0.2065536 0.2168333
## 
## $Covariation_components
##             Total         A         C         E      A(%)
## X <-> X 1.1199492 0.2900561 0.4113080 0.4185850 0.2589905
## X <-> Y 0.7695746 0.3712059 0.1675600 0.2308088 0.4823520
## Y <-> Y 1.0069377 0.5806135 0.2079866 0.2183376 0.5766131
##              C(%)      E(%)
## X <-> X 0.3672559 0.3737536
## X <-> Y 0.2177306 0.2999174
## Y <-> Y 0.2065536 0.2168333
## 
## $All_correlations
##               r_A       r_C       r_E     Total
## X <-> X 1.0000000 1.0000000 1.0000000 1.0000000
## X <-> Y 0.9045451 0.5728871 0.7634776 0.7246867
## Y <-> Y 1.0000000 1.0000000 1.0000000 1.0000000

Output - model fit

# Model fit
model2 <- ref_models(model2)

fit_stats(model2)
##   model      base ep minus2LL  df      AIC       BIC
## 1   ACE Saturated 11 1793.504 789 215.5039 -2386.869
##         CFI       TLI      RMSEA   diffLL diffdf        p
## 1 0.9948431 0.9963598 0.02814297 19.69289     17 0.290238

Vignette

OpenMx

A model (mxModel):

  • matrices (mxMatrix)
  • matrix algebra (mxAlgebra)
  • data (mxData)
  • fit function (mxFitFunction…)
  • expectation (mxExpectation…)
  • sub-models (mxModel)
  • CIs (mxCI), constraints (mxConstraint)

Univatiate ACE

# Twin models in OpenMx
library(OpenMx)

# Univariate ACE model
model <- mxModel(
  name = 'ACE',
  
  # Means
  mxMatrix(type = 'Full',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'mean'),
  mxAlgebra(cbind(mean, mean),
            name = 'expMeans'),
  
  # Variance components
  mxMatrix(type = 'Lower',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'a'),
  mxMatrix(type = 'Lower',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'c'),
  mxMatrix(type = 'Lower',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'e'),
  mxAlgebra(t(a) %*% a, name = 'A'),
  mxAlgebra(t(c) %*% c, name = 'C'),
  mxAlgebra(t(e) %*% e, name = 'E'),
  
  # Expected covariance
  mxAlgebra(rbind(cbind(A + C + E, A + C    ),
                  cbind(A + C,     A + C + E)),
            name = 'expCovMZ'),
  mxAlgebra(rbind(cbind(A + C + E, .5%x%A + C),
                  cbind(.5%x%A + C, A + C + E)),
            name = 'expCovDZ'),
  
  # MZ submodel
  mxModel(name = 'MZ',
          mxData(observed = mz_data,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovMZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'X2')),
          mxFitFunctionML()),
  
  # DZ submodel
  mxModel(name = 'DZ',
          mxData(observed = dz_data,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovDZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'X2')),
          mxFitFunctionML()),
  
  mxFitFunctionMultigroup(c('MZ','DZ')),
  
  # Output tables
  mxAlgebra(A + C + E, name = 'V'),
  mxAlgebra(cbind(V, A, C, E),
            dimnames = list('X', c('Total', 'A', 'C', 'E')),
            name = 'Raw_variance'),
  mxAlgebra(cbind(A / V, C / V, E / V),
            dimnames = list('X', c('A(%)', 'C(%)', 'E(%)')),
            name = 'Variance_components')
)

model <- mxRun(model)

summary(model)
## Summary of ACE 
##  
## free parameters:
##            name matrix row col    Estimate  Std.Error A
## 1 ACE.mean[1,1]   mean   1   1 -0.01540773 0.06578840  
## 2    ACE.a[1,1]      a   1   1  0.73470965 0.08905756  
## 3    ACE.c[1,1]      c   1   1  0.56030417 0.12403417  
## 4    ACE.e[1,1]      e   1   1  0.41025897 0.02890134  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:              4                    396
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:               984.723
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/400
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:        192.723                992.723
## BIC:      -1113.411               1005.916
##       |  Sample-Size Adjusted
## AIC:                 992.9281
## BIC:                 993.2438
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:50:49 
## Wall clock time: 0.04244709 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Univariate output

# Univariate output
mxEval(Raw_variance, model)
##      Total         A         C         E
## X 1.022051 0.5397983 0.3139408 0.1683124
mxEval(Variance_components, model)
##        A(%)      C(%)     E(%)
## X 0.5281518 0.3071673 0.164681

Bivariate ACE

# Bivariate ACE model
model2 <- mxModel(
  name = 'ACE',
  
  # Means
  mxMatrix(type = 'Full',
           nrow = 1, ncol = 2,
           free = TRUE,
           name = 'mean'),
  mxAlgebra(cbind(mean, mean),
            name = 'expMeans'),
  
  # Variance components
  mxMatrix(type = 'Lower',
           nrow = 2, ncol = 2,
           free = TRUE,
           name = 'a'),
  mxMatrix(type = 'Lower',
           nrow = 2, ncol = 2,
           free = TRUE,
           name = 'c'),
  mxMatrix(type = 'Lower',
           nrow = 2, ncol = 2,
           free = TRUE,
           name = 'e'),
  mxAlgebra(t(a) %*% a, name = 'A'),
  mxAlgebra(t(c) %*% c, name = 'C'),
  mxAlgebra(t(e) %*% e, name = 'E'),
  
  # Covariance
  mxAlgebra(rbind(cbind(A + C + E, A + C    ),
                  cbind(A + C,     A + C + E)),
            name = 'expCovMZ'),
  mxAlgebra(rbind(cbind(A + C + E, .5%x%A + C),
                  cbind(.5%x%A + C, A + C + E)),
            name = 'expCovDZ'),
  
  # MZ submodel
  mxModel(name = 'MZ',
          mxData(observed = mz_data2,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovMZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'Y1', 'X2', 'Y2')),
          mxFitFunctionML()),
  
  # DZ submodel
  mxModel(name = 'DZ',
          mxData(observed = dz_data2,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovDZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'Y1', 'X2', 'Y2')),
          mxFitFunctionML()),
  
  mxFitFunctionMultigroup(c('MZ','DZ')),
  
  # Output tables
  mxAlgebra(A + C + E, name = 'V'),
  mxAlgebra(abs(A) + abs(C) + abs(E), name = 'Vabs'),
  mxAlgebra(abs(A) / Vabs, name = 'Aperc'),
  mxAlgebra(abs(C) / Vabs, name = 'Cperc'),
  mxAlgebra(abs(E) / Vabs, name = 'Eperc'),
  mxAlgebra(cbind(diag2vec(Aperc),
                  diag2vec(Cperc),
                  diag2vec(Eperc)),
            dimnames = list(c('X', 'Y'), c('A(%)', 'C(%)', 'E(%)')),
            name = 'Variance_components'),
  
  mxAlgebra(cbind(vech(V),
                  vech(A), vech(C), vech(E),
                  vech(Aperc), vech(Cperc), vech(Eperc)),
            dimnames = list(c('X <-> X', 'X <-> Y', 'Y <-> Y'),
                            c('Total', 'A', 'C', 'E', 'A(%)', 'C(%)', 'E(%)')),
            name = 'Covariation_components'),
  
  mxMatrix(type = "Iden",
           nrow = 2, ncol = 2,
           name = "I"),
  mxAlgebra(sqrt(I * V), name = "SD_total"),
  mxAlgebra(sqrt(I * A), name = 'SD_A'),
  mxAlgebra(sqrt(I * C), name = 'SD_C'),
  mxAlgebra(sqrt(I * E), name = 'SD_E'),
  
  mxAlgebra(solve(SD_total) %&% V, name='r_total'),
  mxAlgebra(solve(SD_A) %&% A, name = 'r_A'),
  mxAlgebra(solve(SD_C) %&% C, name = 'r_C'),
  mxAlgebra(solve(SD_E) %&% E, name = 'r_E'),
  
  mxAlgebra(cbind(vech(r_A), vech(r_C), vech(r_E),
                  vech(r_total)),
            dimnames = list(c('X <-> X', 'X <-> Y', 'Y <-> Y'),
                            c('r_A', 'r_C', 'r_E', 'Total')),
            name='All_correlations')
)

model2 <- mxRun(model2)

summary(model2)
## Summary of ACE 
##  
## free parameters:
##             name matrix row col     Estimate  Std.Error A
## 1  ACE.mean[1,1]   mean   1   1  0.028582935 0.06602392  
## 2  ACE.mean[1,2]   mean   1   2 -0.006697959 0.06398142  
## 3     ACE.a[1,1]      a   1   1  0.229626300 0.19038682  
## 4     ACE.a[2,1]      a   2   1  0.487155814 0.14554958  
## 5     ACE.a[2,2]      a   2   2  0.761981284 0.09837065  
## 6     ACE.c[1,1]      c   1   1  0.525659559 0.07115918  
## 7     ACE.c[2,1]      c   2   1  0.367418707 0.20941586  
## 8     ACE.c[2,2]      c   2   2  0.456054195 0.15874248  
## 9     ACE.e[1,1]      e   1   1  0.417843143 0.02948765  
## 10    ACE.e[2,1]      e   2   1  0.493956895 0.05375951  
## 11    ACE.e[2,2]      e   2   2  0.467266025 0.03298568  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:             11                    789
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:              1793.504
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/800
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:       215.5039               1815.504
## BIC:     -2386.8685               1851.785
##       |  Sample-Size Adjusted
## AIC:                 1816.908
## BIC:                 1816.936
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:50:49 
## Wall clock time: 0.06305766 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Bivariate output

# Bivariate output
mxEval(Variance_components, model2)
##        A(%)      C(%)      E(%)
## X 0.2589839 0.3672615 0.3737545
## Y 0.5766147 0.2065523 0.2168331
mxEval(All_correlations, model2)
##               r_A      r_C      r_E     Total
## X <-> X 1.0000000 1.000000 1.000000 1.0000000
## X <-> Y 0.9045492 0.572894 0.763478 0.7246873
## Y <-> Y 1.0000000 1.000000 1.000000 1.0000000

Bivariate output

mxEval(Covariation_components, model2)
##             Total         A         C         E      A(%)
## X <-> X 1.1199498 0.2900490 0.4113145 0.4185863 0.2589839
## X <-> Y 0.7695757 0.3712036 0.1675628 0.2308093 0.4823484
## Y <-> Y 1.0069384 0.5806155 0.2079854 0.2183375 0.5766147
##              C(%)      E(%)
## X <-> X 0.3672615 0.3737545
## X <-> Y 0.2177341 0.2999176
## Y <-> Y 0.2065523 0.2168331

mlth.data.frame

library(mlth.data.frame)

# mlth.data.frame - multiheaded data.frame
model2 <- multivariate_ace(
  data = data2,
  zyg = 'zyg',
  ph = c('X', 'Y')
)
model2 <- def_ci(model2, ci = 'Variance_components')
model2 <- mxRun(model2, intervals = TRUE)

get_output_tables(model2)
## $Variance_components
##   A(%)-------------------------- C(%)---------------------------- E(%)-------------------------
##   est       lCI        uCI       est       lCI          uCI       est       lCI       uCI      
## X 0.2589905 0.02667819 0.6127657 0.3672559 4.764229e-02 0.5861756 0.3737536 0.2791622 0.4919170
## Y 0.5766131 0.31274212 0.8272015 0.2065536 5.511957e-33 0.4481345 0.2168333 0.1594465 0.2973594
## 
## $Covariation_components
##             Total         A         C         E      A(%)
## X <-> X 1.1199492 0.2900561 0.4113080 0.4185850 0.2589905
## X <-> Y 0.7695746 0.3712059 0.1675600 0.2308088 0.4823520
## Y <-> Y 1.0069377 0.5806135 0.2079866 0.2183376 0.5766131
##              C(%)      E(%)
## X <-> X 0.3672559 0.3737536
## X <-> Y 0.2177306 0.2999174
## Y <-> Y 0.2065536 0.2168333
## 
## $All_correlations
##               r_A       r_C       r_E     Total
## X <-> X 1.0000000 1.0000000 1.0000000 1.0000000
## X <-> Y 0.9045451 0.5728871 0.7634776 0.7246867
## Y <-> Y 1.0000000 1.0000000 1.0000000 1.0000000

M <- get_output_tables(model2)$Variance_components

M[['A(%)']]
##   est       lCI        uCI      
## X 0.2589905 0.02667819 0.6127657
## Y 0.5766131 0.31274212 0.8272015

# Formatting the table for html or latex report
# Powered by knitr and kableExtra
library(kableExtra)

M %>% 
  behead() %>%
  kable(digits = 3,
        align = 'c') %>%
  add_complex_header_above(M)
A(%)
C(%)
E(%)
est lCI uCI est lCI uCI est lCI uCI
X 0.259 0.027 0.613 0.367 0.048 0.586 0.374 0.279 0.492
Y 0.577 0.313 0.827 0.207 0.000 0.448 0.217 0.159 0.297 

# Saving the table for the output
M %>%
  register_output(
    name = 'Table 1',
    caption = 'Table 1: Variance components')
##   A(%)-------------------------- C(%)---------------------------- E(%)-------------------------
##   est       lCI        uCI       est       lCI          uCI       est       lCI       uCI      
## X 0.2589905 0.02667819 0.6127657 0.3672559 4.764229e-02 0.5861756 0.3737536 0.2791622 0.4919170
## Y 0.5766131 0.31274212 0.8272015 0.2065536 5.511957e-33 0.4481345 0.2168333 0.1594465 0.2973594

# Run at the end of the script
# to write all output tables into a single file
# openxlsx package required
write.xlsx.output(file = 'output.xlsx')

Vignettes